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Abstract 

The stability of difference schemes for, in general, hyperbohc systems 
of conservation laws with source terms are studied. The basic approach 
is to investigate the stability of a non-linear scheme in terms of its cor- 
responding scheme in variations. Such an approach leads to application 
of the stability theory for linear equation systems to establish stability 
of the corresponding non-linear scheme. It is established the notion that 
a non-linear scheme is stable if and only if the corresponding scheme in 
variations is stable. 

A new modification of the central Lax-Friedrichs (LxF) scheme is de- 
veloped to be of the second order accuracy. A monotone piecewise cubic 
interpolation is used in the central schemes to give an accurate approxi- 
mation for the model in question. The stability of the modified scheme are 
investigated. Some versions of the modified scheme are tested on several 
conservation laws, and the scheme is found to be accurate and robust. 

As applied to hyperbolic conservation laws with, in general, stiff source 
terms, it is constructed a second order nonstaggered central scheme based 
on operator-splitting techniques. 

1 Introduction 

We are mainly concerned with the stability of difference schemes for hyper- 
bolic systems of conservation laws with source terms. Such systems are used 
to describe many physical problems of great practical importance in magneto- 
hydrodynamics, kinetic theory of rarefied gases, linear and nonlinear waves, vis- 
coelasticity, multi-phase flows and phase transitions, shallow waters, etc. (see. 
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e.g., [B], [in], [H], [21], [301, [32]> El, [37], [H], [12]). We will consider a system 
of hyperbolic conservation laws written as follows (e.g., [17], [32] ) 



where x ={xi,X2, ■ ■ ■ ,xjv}"^ G K^, u = {ui,U2, ■ ■ ■ ,UAf}^ is a vector-valued 
function from x [0, +oo) into an open subset fiu C M*^, f,- (u) — {fij (u) , 
/2j (u) , . . . , fMj (u)}"^ is a smooth function (flux- function) from il^ into M*'^, 
q (u) = {qi (u) , g2 (u) , . . . ,qM (u)}^ denotes the source term, r > denotes 
the stiffness parameter, u*' (x) is of compact support. We will assume that 
T = const without loss of generality. In what follows ||M||p denotes the matrix 

norm of a matrix M induced by the vector norm ||v||p = {J^i and ||M|| 

denotes the matrix norm induced by a prescribed vector norm. M denotes the 
field of real numbers. 

For studying stability and monotonicity of non-linear schemes, the well 
known notion of total variation diminishing (TVD, see, e.g., [T7], [32]) turns 
out to be an useful tool. Actually, the following property 



||AA(v + 5v)-A/'(v)|| < {l + aAt)\\Sv\\ (2) 

is sufficient for stability of a two-step method [32,, however it is, in general, 
difficult to obtain. Here At denotes the time increment, a is a constant inde- 
pendent of At as At ^ 0, V and i5v are any two grid functions {5v will often 
be referred to as the variation of the grid function v), Af denotes the scheme 
operator. At the same time, the stability of linearized version of the non-linear 
scheme is generally not suflicient to prove convergence [20], [32]. Instead, the 
TV-stability adopted in [20] (see also [321 s. 8.3.5]) makes it possible to prove 
convergence (to say, TV- convergence) of non-linear scalar schemes with ease. 
However, the TVD property is a purely scalar notion that cannot, in general, 
be extended for non-linear systems of equations, as the true solution itself is 
usually not TVD [17], [32]. Moreover, one can see in "8", pp. 1578-1581] that a 
TVD scheme can be non-convergent in, at least. Loo, in spite that the scheme 
is TV-stable. Such a phenomenon is, in all likelihood, caused by the fact that 
TV is not a norm, but a semi-norm. 

Nowadays, there exists a few methods for stability analysis of some classes 
of nonlinear difference schemes approximating systems of PDEs (see, e.g., |14j . 
[TB] . [32], [33 , [3H], [H] and references therein). It is noted in [TO] that the 
problem of stability analysis is still one of the most burning problems, because 
of the absence of its complete solution. In particular, as noted in [T3] in this 
connection, the vast majority of difference schemes, currently in use, have still 
not been analyzed. LeVeque [32] noted as well that, in general, no numerical 
method for non-linear systems of equations has been proven to be stable. There 
is not even a proof that the first-order Godunov method converges on general 
systems of non-linear conservation laws [3^ p. 340]. Thus, a different approach 
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to testing scheme stability must be adopted to prove convergence of non-linear 
schemes for systems of PDEs. The notion of scheme in variations (or variational 
scheme [S], [5]) has, in all hkelihood, much potential to be an effective tool 
for studying stability of nonlinear schemes. Such an approach goes back to 
the one suggested by Lyapunov (1892), namely, to investigate stability by the 
first approximation. This idea has long been exploited for investigation of the 
stability of motion [15] . An approach to investigate non-linear difference schemes 
for monotonicity in terms of corresponding variational schemes was suggested 
in [5], [H]. The advantage of such an approach is that the variational scheme 
will always be linear and, hence, enables the investigation of the monotonicity 
for nonlinear operators using linear patterns. It is proven for the case of explicit 
schemes that the monotonicity of a variational scheme will guarantee that its 
original scheme will also be monotone [8]. We establish the notion that the 
stability of a scheme in variations is necessary and sufficient for the stability of 
its original scheme (see Section [21 Theorem |4]) . 

An extensive literature is devoted to central schemes, since these schemes 
are attractive for various reasons: no Riemann solvers, characteristic decompo- 
sitions, complicated flux splittings, etc., must be involved in construction of a 
central scheme (see, e.g., [S], [30], [3T], [32], [H], [13] and references therein) , and 
hence such schemes can be implemented as a black-box solvers for general sys- 
tems of conservation laws [30]. Let us, however, note that the numerical domain 
of dependence [321 p. 69] for a central scheme approximating, e.g., a scalar trans- 
port equation coincides with the numerical domain of dependence for a standard 
explicit scheme approximating diffusion equations [321 P- 67]. Such a property 
is inherent to central schemes in contrast to, e.g., the first-order upwind schemes 
[32l p. 73]. Hence, central schemes do not satisfy the long known principle (e.g., 
[21 p. 304]) that derivatives must be correctly treated using type-dependent dif- 
ferences, and hence there is a risk for every central scheme to exhibit spurious 
solutions. The results of simulations in [39] can be seen as an illustration of 
the last assertion. Notice, all versions of the, so called, Nessyahu-Tadmor (NT) 
central scheme, in spite of sufficiently small CFL (Courant-Friedrichs-Lewy [32]) 
number [Cr — 0.475), exhibit spurious oscillations in contrast to the second- 
order upwind scheme {Cr — 0.95). The first order, 0(Ai -I- Aa;), LxF scheme 
exhibits the excessive numerical viscosity. Thus, the central scheme should be 
chosen with great care to reflect the true solution and to avoid significant but 
spurious peculiarities in numerical solutions. 

Let us note that LxF scheme - the forerunner for central schemes [5], [30] 
- does not produce spurious oscillations. While, from the pioneering works of 
Nessyahu and Tadmor [39] and on, the higher order versions of LxF scheme can 
produce spurious oscillations. The reason has to do with a negative numerical 
viscosity introduced to obtain a higher order accurate scheme (for more details, 
see Section U]). Let us note that there is a possibility to increase the scheme's 
order of accuracy, up to 0((At)^ -f (Ax)^), by introducing an additional non- 
negative numerical viscosity into the scheme. Such an approach is similar to 
the vanishing viscosity method [T7], [32], and hence possesses its advantages, yet 



3 



it appears to be free of the disadvantages of this method, since the additional 
viscosity term is not artificial. With this approach, the second order scheme is 
developed in Section 21 where sufficient conditions for stability of the scheme 
are found. The scheme is tested on several conservation laws in Section [51 

A stable numerical scheme may yield spurious results when applied to a stiff 
hyperbolic system with relaxation (see, e.g., [T], [Ij, [B], [TU], [H], [H], [H], 
[45]). Specifically, spurious numerical solution phenomena may occur when un- 
derresolved numerical schemes (i.e., insufficient spatial and temporal resolution) 
are used (e.g., [T], [H], [IB], [37]). However, during a computation, the stiff- 
ness parameter may be very small, and, hence, to resolve the small stiffness 
parameter, we need a huge number of time and spatial increments, making the 
computation impractical. Hence, we are interested to solve the system, ([1]), with 
underresolved numerical schemes. It is significant that for relaxation systems a 
numerical scheme must possess a discrete analogy to the continuous asymptotic 
limit, because any scheme violating the correct asymptotic limit leads to spu- 
rious or poor solutions (see, e.g., [TU], [M], [13, [37], [II])- Most methods for 
solving such systems can be described as operator splitting ones, [TT], or meth- 
ods of fractional steps, [6]. After operator splitting, one solves the advection 
homogeneous system, and then the ordinary differential equations associated 
with the source terms. As reported in [TB], this approach is well suited for the 
stiff systems. We are mainly concerned with such an approach in Section 

2 Stability of difference schemes 

Let us consider the following non- linear explicit scheme arising, e.g., in numerical 
analysis of nonlinear PDE systems: 

v^+i =Hr(vl\v^\...,v^'), :r!„CR^^R^o, z e ^i, n,n+lGuj2, (3) 

where v" G M.^" denotes a vector- valued grid function, N = NqI, i G cji denotes 
a node of the grid uji = {1,2,...,/}, n & uj2 denotes a node (time level) of the 
grid = {0, 1, . . . , A/}, = {Hl]^, Hf.^, H^^^J^ is a vector-valued func- 
tion with the domain and range belonging to and R^° , respectively. Notice, 
H" depends also on scheme parameters (e.g., space and time increments), how- 
ever, this dependence is usually not included in the notation. We will assume 
that n in ([3]) denotes the time level, i„ (= nAt). Thus, the time increment will 
be represented by At = t^a.x/M — const, where tmax denotes some finite time 
over which we wish to compute. If we introduce the additional notation 

v"- {(vD^ , (v^)^ , . . . , (VlfY , H"= {{U^f , (H^)^ , . . . , {UrfY , (4) 

then the scheme ([3]) can be written in the form 

v"+i=H"(v"), H":aiCR^^R^, n, n+ 1 e = {0, 1, . . . , M} . (5) 

As usual (e.g., [40;, p. 62]), for mappings f : f^/ C R^ ^ R^ and g : C 
R^, the composite mapping h = g o f is defined by h (v) — g (f (v)) for 
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all V gJI/i — {v Gftf I f (v) E fig}. Using the composite mapping approach, we 
rewrite Scheme ([5]) to read 

y = F(x), F : C -> R^, (6) 

where the following notation is used: x = v°, y = v*-'^, F = H^^~^ oH^-'^^^ o . . .o 
HO, ftp = {vOeOo I vi = HO (vO) e r!i I . . . I v^~^ = B.^'-^ {v^'-^) e I^m-i}. 
Let the scheme parameters (including time increments) be represented by a vec- 
tor s belonging to some normed space with the norm |s|. 

Since differentiability of H" as well as F will be a key element in the fol- 
lowing, let us note that the composite mapping F will be Frechet-differentiable 
[101 item 3.1.5] if all of the maps, H", are Frechet-differentiable [IHl item 3.1.7]. 
However, if all of the maps are Frechet-differentiable, but one that Gateaux- 
differentiable 00] item 3.1.1], then the composite mapping F has a Gateaux- 
derivative [40l item 3.1.7]. Notice, if there exist at least two maps having 
Gateaux-derivatives, then F need not be differentiate [40l E 3.1-7]. 

Scheme ^ is said to be stable (see, e.g., |T4], [17], 02], [l^, [48], [49]) if 
there exist positive so, C — const such that for all x, x* G fi^? the following 
inequality is valid 

||F(x,)-F(x)|| <C||x,-x||, Vs: |s| < sq- (7) 

Thus, Scheme (jG]) will be stable ijf (if and only if) the function F will be 
Lipschitz for a constant C. 

To be more specific, let us consider the "slit plane" [22l in polar coordinates 
{r,0) 

np = {(r, 6*) I < r < oo, -tt < S < tt} C M^ (8) 
and the function F = {Fi (r, 9) , F2 (r, 6*)}^ such that [22] 

Fi =r, F2 = e/2. (9) 

If we take {r,9)^ — (ro,— tt + e), {r,0) = (ro,7r — e), and tq = const, then 
obviously the mapping ^ is not Lipschitz, since C in ([7]) tends to infinity as 
e — » 0. Therefore, we have to conclude, in view of the above definition, that 
Scheme © is not stable, even though the function F, ([5]), is locally Lipschitz 
for C = 1, and, further, F is the non-stretching mapping of the "slit plane" 
(|8]) into the right semi-plane. Hence, the preceding definition of stability needs 
minor changes. 

A set C R-'^ is said to be path-connected if every two points x, x* G 
can be joined by a continuous curve (7 : [0, 1] C R ^ fi, [22], [Ml P- 113]) of 
finite length, ^(7). The intrinsic metric [32] An in a path-connected set 17 is 
defined as 

An(x,x*) = inf L(7), 7: x=7(0), X* =7(1), L(7) <cx). (10) 
7CSI 

An open ball (of radius r) about x SR^ is denoted by B (x, r) (or just i?x)- 
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Definition 1 Let Qp in (0) be path- connected. Scheme is said to be stable 
if there exist positive sq, C ~ const such that the following inequality holds 

||F(x,)-F(x)|| <CAn^(x,x,), y^,^,enF, V s : |s| < so- (h) 

Notice, Scheme © is stable, since Inequality (ITT|) holds for C — I. 

Lemma 2 Let the path- connected ilp of (0j be open in R^. Scheme (0j will be 
stable in terms of Definitionl^iS F in (0j will be locally Lipschitz for a common 
constant C, for all scheme parameters s such that |s| < sq. 

Proof. Suppose Scheme ^ is stable, i.e. (|lip is valid. Choose any point 
X GfJiT'. Since fi^? is open, there exists a radius r such that B{x,r) C ftp- 
Choose any point x, G i? (x, r), and let 7^ be the straight line segment joining 
the points x, x^, £ i? (x, r) . In view of pT|) , F in ([S]) will be locally Lipschitz for a 
common constant C, for all s : |s| < sq, since An^ (x, x») = L (7^) = ||x, — x||. 

Conversely, suppose that F in ([S]) is locally Lipschitz for a common constant 
C, for all s : |s| < sq- Let some points x, x, (z flp he joined by a continuous 
curve 7. In view of pO)) . the curve 7 can be taken such that L (7) < Aoj^ (x, x»)+ 
e for an arbitrary e > 0. Given any point z G 7, there is a ball C ftp. The 
balls {Bz} form an open cover of 7. Since the mapping 7 : [0, 1] C K ^ is 
continuous, the curve 7 is compact [29l p. 94]. Hence, by the compactness of 7, 
{B2} has a finite subcover consisting of balls B^ = -Bzi, ^za: • ■ -i ^za- = B^,. 
Since F is locally Lipschitz, we find 

||F(zfc+i)-F(zfc)|| <C7||zfc+i-Zfc||, fc = l,2,...X-l, Vs: |s| < sq. (12) 

Then, by virtue of p^ . we find 



|F(x,)-F(x) 



^[F(z,+i)-F(zfc) 

k 



<C^||zfc+i-Zfc|I < 

k 



CL{j) <C An A^,^*)+£C, Vs:|s|<so- (13) 

By letting e ^ 0, we find that (fTTj) holds. ■ 

Let us find the necessary and sufficient conditions for the stability of Scheme 
©. Let W^-°° (np) denote the Sobolev space, and let F = {Fi , i^2, ■ ■ • , ^V}^ 
in Then, F^, i = 1,2, . . . , N, (and, hence, F) is locally Lipschitz (in the 
sense of having representatives) iff Fi G W^'°° (ftp) (see, e.g., [211 Theorem 
4.1]). Let VFi denote the distributional gradient of Fi, and let SF, Sx G 
denote variations. The following equality 

5F = F' • (5x, F' = {V^^i , \/F2, \7Fn}^ , (14) 

will be viewed as the scheme in variations for (161). 



Lemma 3 Linear Scheme ^4\ ) stable iff there exist positive sq, C = 

const such that 

\\F'\\ < C = const, VxeQp, Vs: js] < sq. (15) 
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Proof. The sufficiency is obvious. Actually, by virtue of (|15p . we find that 
||<5F|| = ||F'.<5x||<|lF'||||<5x||<C||5x||,i.e. 

||(5F|| < C ||(5x|| , Vxel^F, Vs: |s| < sq. (16) 

Conversely, suppose that is valid. Then, in view of [231 Theorem 2, p. 224], 
we write 

l|F' -^xll MFII C||<5x|| ^ 

IIF'II = sup " " = sup ^ < sup = C. (17) 

Hence, ([15]) holds, in view of p7|) ■ 

Theorem 4 Consider Scheme 0). Lef t/ie path- connected be open, F = {i^i , 
i^2, i^Tv}^ &e ftownrferf, and let ||F'|| = Hf^ll, = {||VFi||^ , IIVF2I' 
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. . . , ||Vi^Ar||j^}"^, Vi<i, i = 1, 2, . . . , N, denote the distributional gradient of Fi. 
Then, Scheme (0) will be stable iff its scheme in variations, wi/l be stable. 



Proof. The proof is trivial. Actually, Scheme ^ is stable F is locally 
Lipschitz for a common constant C (Lemma [2]) <;=^ Fi G W^'°° (ilp) (see [2^ 
Theorem 4.1]) (fT5| holds Scheme in variations, (fT4| . is stable (Lemma 
[3]). ■ 

Notice, if F in (l6|) is Gateaux-differentiable, then VF^ (see Theorem |4|) de- 
notes the classical gradient of Fi, and, hence, it may be taken that ip — {|| VFi | 
, ||Vi^2|| , ■ • ■ , ||Vi^Ar||}^, see also [3 Theorem 3]. 

3 Monotone piecewise cubics in construction 
of central schemes 

In this section we consider some theoretical aspects for high-order interpola- 
tion and employment of monotone piecewise cubics (e.g., [T2], [23) in con- 
struction of monotone central schemes. We will consider explicit schemes on 
a uniform grid with time step At and spatial mesh size Aa:, as applied to the 
following hyperbolic 1-D equation 

^ + — f(u) =0, <„ <t<i„+i =i„ + Ai, u(x,t„) = u"(a;), (18) 
Using the central differencing, we write 



du 

'di 



ri+0.5 _ n 



fri+0.25 fn+0.25 

dx 



h+1 h 



{{Axf) . (20) 
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By virtue of P^ - ipO]) we approximate ([T5|) on the cell [a;^, x^+i] x [in, in+o.s] by 
the following difference equation 

^i+0.5 — ^1+0.5 Si j • i^J-; 

As usual, the mathematical treatment for the second step (i.e., on the cell 
[xi_o.5, Xi+o.s] X [tn+o.5,tn+i]) of a Staggered scheme will, in general, not be 
included in the text, because it is quite similar to the one for the first step. 

Considering that (PT|) approximates (|18p with the accuracy 0((Aa;)^+(At)^), 
the next problem is to approximate vf_)_g 5 and g""*""'^^ in such a way as to retain 
the accuracy of the approximation. For instance, the following approximations 

vr+o.5-o.5(vr+vr+o + o((Axf), gr°-^5 = f(vn+o(Ai), (22) 

leads to the staggered form of the famed LxF scheme that is of the first-order 
approximation (see, e.g., [171 P- 170]). One way to obtain a higher-order scheme 
is to use a higher order interpolation. At the same time it is required of the 
interpolant to be monotonicity preserving. Notice, the classic cubic spline does 
not possess such a property (see Figure [TJi). Let us consider the problem of 
high-order interpolation of v"_(.q 5 in with closer inspection 

Let p = p (a;) = {p^ (x) , . . . (x) , . . . ,p™ {x)}'^ be a component-wise 
monotone piecewise cubic interpolant (e.g., [T^], [21])) and let 



P^^P {X^) , P- = p' (Xi) , Ap, = p,+i - p,. 



(23) 



where p[ denotes the derivative of the interpolant at x = Xi. The diagonal 
matrices and in are defined as follows 

A, = dtag {alal <} , 1, = dtag {f]], (3l . . . , /Sf} ■ (24) 

The cubic interpolant, p = p{x), is component-wise monotone on [xi,a;i+i] iff 
one of the following conditions (e.g., [H], [55]) is satisfied: 

(af - 1)' + {a'y 1) _ 1) + (;3f - 1)' - 3 («f + - 2) < 0, (25) 

af + /3f < 3, > 0, /3f > 0, Vi, k. (26) 

As reported in [28], the necessary and sufficient conditions for monotonicity of 
a piecewise cubic interpolant originally given by Ferguson and Miller (1969), 
and independently, by Fritsch and Carlson ^12. . The region of monotonicity 
is shown in Figure [T|d. The results of implementing a monotone piecewise 
cubic interpolation when compared with the classic cubic spline interpolation, 
are depicted in Figure [T^. We note (Figure [T^) that the constructed function 
produces monotone interpolation and this function coincides with the classic 
cubic spline at some sections where the classic cubic spline is monotone. 
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(a) (b) 




Figure 1: Monotone piecewise cubic interpolation, (a) Interpolation of a 1-D 
tabulated function. Circles: prescribed tabulated values; Dashed line: classic 
cubic spline; Solid line: monotone piecewise cubic, (b) Necessary and suffi- 
cient conditions for monotonicity. Horizontal hatching: region of monotonicity; 
Unshaded: cubic is non-monotone. 



Using the cubic segment of the piecewise cubic interpolant, p — p(a;), 
(see, e.g., [12], [28]) for x G [xi,Xi+i], we obtain the following interpolation 
formula 

p,+o.5 = 0.5 (p, + p,+i) - ^ (p^+i -p'i}+0 {{AxY) . (27) 

If p (x) has a continuous fourth derivative, then r = 4 in ((7f)) . see e.g. [IT] p. 
111]. However, the exact value of p^ in ([?/)) is, in general, unknown, and hence 
to construct numerical schemes, employing formulae similar to (j27p . the value 
of derivatives p^ must be estimated. 

Using ([27]) and the second formula in ([22]) we obtain from (I2ip the following 
scheme 



-Uo'i = 0.5 K + vr+0 - ^ (div, - dD - ^^^^H/^, (28) 

where d" denotes the derivative of the interpolant at a; = Xi. In view of (|27p 
and the second formula in (|22p . the local truncation error [32J p. 142], ip, on a 
sufficiently smooth solution u(x, t) to is found to be 

i; = (At) + O (^^) + O ((At)2 + (Axf) . (29) 

In view of (HHj) we conclude that the scheme ([^5]) generates a conditional approx- 
imation, because it approximates only if {AxY /"At — > as Ax and 
At — > 0. Let d" be approximated with the accuracy O ((Aa;)^), then the value 
of r in (|29p can be calculated (see Section [6l Proposition [5|) by the following 
formula 

r = min(4, s+ 1). (30) 
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Interestingly, since (^5]) provides the conditional approximation, the order of 
accuracy depends on the pathway taken by Ax and At as Ax — > and At 0. 
Actually, there exists a pathway such that At is proportional to (Ax)^ and 
the CFL condition is fulfilled provided /i > 1 and Ax < Axq, where Axq is 
a positive value. If we take n = 1 and s > 1, then we obtain from that 
the scheme (1^51) is of the first-order, li ji — 2 and s > 3, then is of the 
second-order. However, if /i = 2 and s = 2, then, in view of ((29| and (|30p. the 
scheme (pS)) is of the first-order. Moreover, under fi — 2 and s — 2, the scheme 
will be of the first-order even if g""*"*^'^^ in (PT|) will be approximated with the 
accuracy 0{{AtY). It seems likely that Example 6 in [30] can be seen as an 
illustration of the last assertion. The Nessyahu-Tadmor (NT) scheme with the 
second-order approximation of d" is used [30] to solve a Burgers-type equation. 
Since At — 0{{Ax) ) [30 , the NT scheme is of the first-order, and hence it 
can be the main reason for the scheme to exhibit the smeared discontinuity 
computed in [30l Fig. 6.22]. 

The approximation of derivatives p.^ can be done by the following three 
steps fT2|: (i) an initialization of the derivatives p^; (ii) the choice of subregion 
of monotonicity; (iii) modification of the initialized derivatives p^ to produce a 
monotone interpolant. 

The matter of initialization of the derivatives is the most subtle issue of this 
algorithm. Actually, the approximation of p[ must, in general, be done with 
accuracy 0{{Ax)^) to obtain the second-order scheme when At is proportional 
to (Ax) , inasmuch as central schemes generate a conditional approximation. 
Thus, using the two-point or the three-point (centered) difference formula (e.g. 
[28] . [41]) we obtain, in general, the first-order scheme. The so called limiter 
functions [23 lead, in general, to a low-order scheme as these limiters are often 
0{Ax) or 0((Aa;)^) accurate. Performing the initialization of the derivatives p[ 
in the interpolation formula P7)) by the classic cubic spline interpolation |i5] . 
we obtain the approximation, which is 0{{Ax)^) accurate (e.g., [17], [IH]), and 
hence, in general, the second-order scheme. The same accuracy, 0{{Ax)^), can 
be achieved by using the four-point approximation |28) . However, the efficiency 
of the algorithm based on the classic cubic spline interpolation is comparable 
with the one based on the four-point approximation, as the number of multi- 
plications and divisions (as well as additions and subtractions) per one node 
is approximately the same for both algorithms. We will use the classic cubic 
spline interpolation for the initialization of the derivatives P'j in the interpola- 
tion formula l|27p. as it is based on the tridiagonal algorithm, which is 'the rare 
case of an algorithm that, in practice, is more robust than theory says it should 
be' [ig. 

Obviously, for each interval [xi, Xi+i] in which the initialized derivatives p,^, 
p-_i_i such that at least one point (af, /3^) does not belong to the region of 
monotonicity ([^ - ([^ . the derivatives p-, p'i_^_i must be modified to p-, p-_,_]^ 

such that the point (3, , /3j ) will be in the region of monotonicity. The modifi- 
cation of the initialized derivatives, would be much simplified if we take a square 
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as a subregion of monotonicity. In connection with this, we will make use the 
subregions of monotonicity represented in the following form: 

< a,^ < 4N, 0</3'^ < 4K, k, (31) 



where H is a monotonicity parameter. Obviously, the condition (|3ip is sufficient 
for the monotonicity (see Figure [T]d) provided that < H < 0.75. 

Let us now find necessary and sufficient conditions for (|27p to be monotonic- 
ity preserving. By virtue of (j23p , the interpolation formula (j27p can be rewritten 
to read 



p,+o.5 = [0-5I + g j ■ p^ + (^0.51 - g ' j • p^+i. (32) 

The coefficients of (|32l) will be non-negative iff — ai\ < 4. Hence (|27| will 
be monotonicity preserving iff ([31]) will be valid provided < H < 1. Notice, 
there is no any contradiction between the sufficient conditions, (|5T|) provided 
< K < 0.75, for the interpolant, p = p(a;), to be monotone through the 
interval [xi,Xi+i], and the necessary and sufficient conditions, ([3T|) provided 
< H < 1, for the scheme (15^ to be monotonicity preserving. In the latter case 
the interpolant, p = p(a;), may, in general, be non- monotone, however at the 
point i + 0.5 the value of an arbitrary component of Pi+o.s will be between the 
corresponding components of Pi and Pi+i. 

To fulfill the conditions of monotonicity (pT|) , the modification of derivatives 
p^ = , . . . can be done by the following algorithm suggested, in 

fact, by Fritsch and Carlson [T^ (see also P5]): 

5f 4Nminmod(A*^_i, Af), := minmod(pf , Sf ), N = const, (33) 

where A*^ — {p^^i — pf)/Aa;, the function minmod(a;, y) is defined (e.g., [55], 
[3D], [35], [H], [OTI) as follows 

minmod(x,?/) = — [sgn{x) + sgn{y)]min{\x\,\y\). (34) 

Let us note that instead of point values, v"_(.q 5, employed in the construction 
of the scheme (pij) . it can be used the cell averages (e.g., [S], [30], [52]) calculated 
on the basis of the monotone piecewise cubics. In such a case we obtain, 
instead of (P7|l . the following interpolation formula 

Ax 

Pi+0.5 = 0.5 {pi + pi+i) - X— (p-+i - P-) , (35) 
where >c — 2/3. The region of monotonicity in this case will also be 

< A, < 4NI, < < 4KI, < H < 1, Vi. (36) 



Notice, the interpolation formula ((35)) coincides with ((27)) under >c = 1. Thus, 
in view of the interpolation formula p5p . the staggered scheme ()21l) is written 
to read 



."+0-5-05(v" +v")-><— (d" d"^ A^ f(v»^0-f(vr) 
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where d" denotes the derivative of the interpolant ai x — Xi, the range of values 
for the parameter x is the segment < x < 1. If x = 1 (or x = 0), then Scheme 
(P7|) coincides with the scheme (1^51) (or with the LxF scheme, respectively) . As 
it was shown above, the scheme ([57)) is of the first order provided At = O (Ace). 
The central scheme (|57)) . approximating the 1-D equation with the first 
order, will be abbreviated to as COSl. 



4 Construction of central schemes 

We will consider explicit schemes on a uniform grid with time step At and spatial 
mesh size Ax. In view of the CFL condition [32], we assume for the explicit 
schemes, that At = O (Ax). Moreover, we will also assume that Aa; = O (At), 
since a central scheme generates a conditional approximation to Eq. ([THj) (see 
Section[3|). In such a case, the following inequalities will be valid, for sufficiently 
small At and Aa;, 

voAt < Ax < HoAt, J^OjMo = const, < < /Xq. (38) 

Notice, for hyperbolic problems it is often assumed that At and Ax are related 
in a fixed manner (e.g., |32l p. 140], |47l p. 120]), i.e. it is assumed that At 
and Aa; fulfill a more strong condition than ([38]) . 

Scheme ([^5]) is of the first-order, 0{At + (Ax)^), and non-oscillatory LxF 
scheme is of the first-order, 0{At + Ax). Let us demonstrate that ([^5]) is, in 
fact, LxF scheme with a negative numerical viscosity added to obtain a higher 
order approximation to Eq. p8|) with respect to x. We rewrite Scheme ([37| to 

''""^ v^^+f-^-vIVo.s , f(vIVi)-f(vr) _ 



O.bAt 

Aa;2 v" - 2v: 



At Aa;2 " AAt Ax ' ^^^^ 

Notice, the second term in the right-hand side of ([55[) is, in fact, the negative 
numerical viscosity. Without this term {k = 0), Scheme ([5^ would be LxF 
scheme. As it is demonstrated in Section [5l Scheme ([25[) can exhibit spurious 
oscillations in contrast to LxF scheme. Interestingly, there is a possibility to 
improve Scheme p7l) by introducing an additional positive numerical viscosity 
such that the scheme's order of accuracy would increase up to 0((Ai)^ -I- (Aa;) ). 
Let us approximate v"^o.5 ^^'^ g,""*""'^^^ in ([^ with the accuracy 0((Aa;) -|- 
(Ai)^). Using Taylor series expansion, we write 



f(vr) 



dt 

By virtue of the PDE system, p8|) . we find 



At 
T 



O (At^) . (40) 



di _ di du _ di di du _ fdiy du 
dt du dt du du dx \ du I dx 



(41) 
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Using the interpolation formula (|35p and the formulae (|^01) - (|¥T|) . we obtain from 
(|2T|) the following second order central scheme 



n+0.5 _ 



0.5 (vIVi 



- X— (dr+i 



Atf(vf+i) -f(vr) 



Ax 



(Ml 

8Aa; 



(A^O'-diVi-(Ar)'-dr 



(42) 



where ^ is introduced by analogy with m: in ([57|) , and hence < ^ < 1 . Scheme 
((42|) coincides with ((37)) provided that ^ = 0. Since df is the derivative of 
the interpolant at x = x^, the last term in right-hand side of (j42p can be seen 
as the non-negative numerical viscosity introduced into the first order scheme 
([37|. Owing to this term, Scheme (j42|) is 0{{Ax)'^ + (At)^) accurate, provided 
that ^ = 1. Thus, we are dealing with the vanishing viscosity method [T7j, [32] 
and, hence, in view of [17, Theorem 3.3], the scheme, (l42l) . satisfies the entropy 
condition. The central scheme (|42)) . approximating the 1-D equation ([18]) with 
the second order, will be abbreviated to as C0S2. 



4.1 Stability of the second-order scheme COS2 

In view of TheoremHJ the stability of (|42]) will be investigated on the basis of its 
variational scheme. It is assumed that the bounded operator A (= di (u) /du) 
in ([T5|) is Frcchet-differentiable on the set Vl^ C M^^, and its derivative is 
bounded on fJu. Hence, the following inequalities are valid 



sup ||A|1 < A„ 



< oo. 



\\6A7\\ = 



aA^ 



9vr 



< a A UK 



(43) 



where Amax, "a 
we write 



const. Considering that v" in (I42|) is Lipschitz-continuous 
< C^Ax, = const. 



11*1 ~ "j+ill — ^v-^-^i — (44) 

By virtue of ([23]) . the second term in right-hand side of ([42]) can be written in 
the form 

(dr+i - d,") = I (Br - K) ■ (vr+i - vr) . (45) 

Then, the variational scheme corresponding to (jl^ is the following 



K+'o°5'=0.5(K + <5vr+i) + - 



+-Dr-(<5vr-<5v^Vi) 



8Ax2 



{ [<5 ((A^+i)' • B^] • HVi - - K^^^')' • "^01 • (^^+1 ^ } + 



'8Aa;2 



(A-Vi) 

At 
2Aa; 



'.B,-(An'-A,] .(<5vr+i-<5vr) + 
(Ar • <5vr - A^+i • K+i) , 



(46) 
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where = diag {D^^, D^2> ■ • ■ > ^"m} = K ~A.f. By virtue of ([Ml), we find 
that -4HI < < mi, and hence' -SHI < SBf < 8KI. Thus, we may write 
that 

pB^II < 8H. (47) 

By virtue of ([25]), (gT]) and (jH]), and since < H < 1, we find the following 
estimation for the second term in right-hand side of (|46p : 



< 



(48) 



By virtue of jMl), dill), (041), and since < ^ < 1 and the CFL number 

Cr = AiAmax/ < 1, wc find the following estimation for the fourth and fifth 
terms in right-hand side of (|15|) : 



•8Ax2 



-5 Ar+i 



(v^Vi - vr) 



</3^At||,Svr||, Pa ^ const, (49) 



where depends on the other constants, namely, on a a, Amax, Cy, (Iq. 

In view of (|48|) and (|49)) . Scheme ((46)) will be stable if the following scheme 
be stable (see [48, pp. 390-392], 7, Theorem 7]). 

6v-+,%' ^ 0.5 (Jvr + <5v^Vi) + >r • {S< - S^7+i) + 



'8Ax2 



At 



We rewrite ([50)) to read 



(50) 



'5vIV+o°5' = 0.5 (I + ED • 5vr + 0.5 (I - Ef+O ■ ^v^+i. 



(51) 



where 



E" 



4 4(Aa;)^ 



(AD' ■ A, 



At 
Ax 



A", 



(52) 



+ 4 ^4(Aa;)' 



(ADi)'-B,-(Ar 



At 

A:c^^+i- 



(53) 



Since the operator A (= 9f (u) /du) is Frechet-differentiable, and its derivative 
is bounded, see (|i5|) . we get, by virtue of ([33]) and [301 Corollary 3.2.4], that 

||EDi - E7II = ^ ||ADi - Aril < l^a^ ||<5vr|| < a^C.At. (54) 
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We find, in view of the first inequality in and (I^Bl) . that the spectrum 
s(Ef) C [-Xe,Xe], where 



f AtX„ 



—. — A,. 



(55) 



V Ax J ' ' Ax' 
Hence, by virtue of jjj, Theorem 7] we find that the scheme ([STj) will be stable if 
max 0.5(|1 + A| + |1-A|) ^ 1, Vi,n. (56) 

We obtain from (|56)) the following condition for the stability of the variational 
scheme (|46)) 

(>^-^a2)H + a < 1, a = ^^^<l. (57) 

Ax 

Thus, in view of Theorem 2] (see also |71 Theorem 3]), the scheme will be 
stable if be valid. 

Let us note that the parameters >^ and ^ are taken as constant in Scheme 
P^ . However, in practice, it can be convenient to take that = >^(v") and 
^" = )■ In such a case the condition, ([ST)) . for the stability of (|42l) can be 
grounded in perfect analogy to the above, if 



ll'5< 



(9> 



where a^, — const. 



<a.ii<5vrii, \\sa\\ = 



< 



«5l|<5vr||, (58) 



4.2 Operator splitting schemes 

By virtue of the operator-splitting idea [5], [TT], [TH], [32] (see also LOS in [iH]). 
the following chain of equations corresponds to the problem ([T]) 

^^ = ^q(U), i„<f<t„+o.5, U(x,t„) =U"(x), (59) 
1 It t 'I 

-J2 9^^^ (U) = 0, i„+o.5 < t < U (x, i„+o.5) = U"+°-' (x) , (60) 



2 at 



where U" (x) denotes the solution to at i = t„, U"+°-^ (x) denotes the 
solution to (|59p at i = t„+o.5. If a high- resolution method is used directly for 
the homogeneous conservation law (|60p . then it is natural to use a high-order 
scheme for (l59l) . As applied to, in general, stiff (r -C 1) System ([T]), the second 
order schemes can be constructed on the basis of operator-splitting techniques 
with ease if (|59l) will be approximated by an implicit scheme and (j60p by an 
explicit one, see Proposition [H] in Section [HI As an example, let us develop 
a central scheme for a 1-D version of ([1]). After operator-splitting, the 1-D 
equation can be represented in the form 



2 ar = 7^(u)^ 



tn <t <tn+0.25, U (x, t„ ) = U" (x) , 



(61) 
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1 dXJ 3 

2^ + ^f(U)=0, i„+o.25 < i < <n+o.5, U(x,t„+o.25) = U"+°-25(a;). (62) 

Let us first consider the case when the following first-order implicit scheme be 
used for ([6T|) 

vr+°-25=vr + ^qK+°-25), (63) 

and a central scheme with nonstaggered grid cells will be used for ([S^ . To 
eliminate the staggering in ([57]) . we can define, e.g. [23], the nonstaggered cell- 
average as the average of its two neighboring staggered cell-averages. Then, by 
virtue of (l37ll. we find 



,ri+0.5 r,r / ri+0.25 , r,, Ji+0.25 , , Ji+0.25\ ^^^^ { -,71+0.25 .n+0.25 
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<-n+0.25 _ fr^+0.25^ /r^N 

It is clear that Scheme (|64| approximates ((62|l with the accuracy 0{At+(Ax) ), 
however, in view of Proposition [5] in Section [SI Scheme (|63p - ([Ml) . taken as a 
whole, is of the second order approximation for the 1-D version of ([T|). 

Let us develop another nonstaggered central scheme approximating a 1-D 
version of ([T]) with the accuracy 0((At)^ -I- {Ax)'^) and such that its components 
(after operator splitting) will be of the second order. It can be done on the 
basis of the second order scheme ([63]) . ([64]) with ease. Actually, adding to and 
subtracting from Equation (|126p (see Section [SI Proposition [H]), rewritten for 
tn < t < t„_|_o.5, the same quantity, we obtain (after operator splitting) the 
following scheme, instead of ([63]) . ([64|l . 



V 



-^^^ T 32 \dt^ j , ' 



{At)' fd^Vy^''-'^ At .^„+o.25 ,„+o.25 



2 / o2tt\ n+0.25 

-0.25 _ fn+0.25\ 

32 V 4Ax ^'^+1 

Thus, Scheme (|65p as well as Scheme (|66p are of the second order, and Scheme 
35|) - ([S5)) . taken as a whole, is of the second order as well. 

Using Taylor series expansion, and central differencing, we find 



.^n-l-0.125 _ ^n+0.25 



At rdv 



1 /AA Va^u^ "+°'' 



2 V 8 / \ dt^ 



at 

0{iAtf), (67) 



n+0.125 



vr°-^^ = vr ^ ) +0i{Atf]. (68) 
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We obtain, by virtue of (gT]), that 



where A =9f /dV. Then 



d_ 

dx 



^ A 



dx 



-1 n+0.25 



_ 1 

~ ~Kx 

n+0.25 



dx 



ri+0.25 



dx 



ri+0.25 



( \n+0.2by 



■n+0.25 



Ax 



O 



(70) 



where {A'l+^ff = 0.5 ((A:!+°-25)2 + [A^l+'^-^^f). By virtue of m, iU-lIZOl), 
we rewrite Scheme (l65])-(l66l) to read 



^n+0.125 _ .^n+0.25 



n+0.25 



At 



n+0.125 



.,n+0.25\ 



At 



,n+0.5 



0.25 (v"+°-25 + 2v: 



n+0.25 



ri+0.25\ 



n+0.125 



(71) 
(72) 



■t+i 



^ 16 



n+0.25 

+ 1 



dn4^.25) 



8 (Ax)' 



(Aiv^f )^ . (vjv^"' - vr°"^) - (Ai^f )' • (vr+°-^^ - vl^-^) 



;'fn+0.25 

4Ax ^ 



pn+0.25\ 
I 



(73) 



where ^ is introduced in the third term in the right-hand side of (|73p by analogy 
with Scheme (|42|) . and, hence, 0<^<1. If^ = 0, then ([73|l coincides with 
([M]) . being 0(At + (Ax)^) accurate. If = 1 and ^ = 1, then Scheme (|7T|) -d73 |) 
approximates the 1-D version of ([T]) with the accuracy 0((A<)^ + (Ax)^). 

By analogy with the scheme C0S2, ((42|) . we find the conditions for the 
stability of Scheme ([75]) using its scheme in variations. Scheme ([75]) will be 
stable if 

o.5^c,2 + a < 1, Cr = < 1. 



Ax 



(74) 



Let us note that in practice (e.g., [32], [48]) the operator-splitting techniques 
find a wide range of application in designing economical schemes for Eq. (|60p 
in the domains of complicated geometry. The resulting method, in general, will 
be only first-order accurate in time because of the splitting [JS]. Thus, in 
line with established practice we will replace the multidimensional Eq. ([5U|) by 
the chain of the one-dimensional equations: 

1 dXJj d 

2N~Qf + Q^^^ (^j) ^ ^' tn+0.5+(j-l)/{2N) < t < tn+O.b+j / {2N) , (75) 
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where (a;, i„+o.5+(j-i)/(2JV)) = U^-i (x, i„+o.5+(j-i)/(2Af)) , j = 1,2, . . .iV, 
Uo(x,t„+o.5) denotes the solution to ((59)) at t — tn+o.5- Eq. ([59]) will be 
approximated by a first-order imphcit scheme or a second-order imphcit Runge- 
Kutta scheme. In particular, it will be used the following Runge-Kutta scheme 

^„+0.25 ^ ^„+0.5 _ (^«+0.5) ^ ^ ^„ ^ (vr+"-'^) , (76) 

since this scheme possesses a discrete analogy to the continuous asymptotic 
limit. 



5 Exemplification and discussion 

In this section, we are mainly concerned with verification of the second order 
central scheme C0S2, (|42]). 



5.1 Scalar non-linear equation 

As the first stage in the verification, we will focus on the following scalar 1-D 
version of the problem ([T|): 

du Q 

^ + ^/(w)=0, a; e M, < t < Tn^ax; u{x,t\^^ = u"" {x) . (77) 

We will solve the inviscid Burgers equation (i.e. / {u) = v?'/2) with the follow- 
ing initial condition 

u {x, 0) = I ^q"' I ^ |J^^^ , hn > h^, no = const + 0. (78) 



The exact solution to ([77|l . (|78|1 is given by 



where T = 2S'//uo, S = hj^ — h^, 



f^uo, hL < X <b, b = UQt + Hl 
ui{x,t) = -^ Uq, b < X < O.duot + hn , (80) 

0, X < Hl or X > 0.5uot + Hr 

{2S{x-hL) u ^ ^ T 

x<hr orx>L ' ^ ^ 

L = 2y/S^ + 0.5uoS {t-T) + hL. (82) 



First, it will be used Scheme (|42l) under ^ = 0, i.e. the first order in time 
central scheme COSl, (j37p . The numerical solutions were computed on a uni- 
form grid with spatial increments of Aa; = 0.01, the velocity uo = 1 in ([78]) . 
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Figure 2: Inviscid Burgers equation. The scheme COSl {>c = 1) versus the 
analytical solution. Crosses: numerical solution; Solid line: analytical solution 
and initial data. C,. = H = 0.5, Aa; = 0.01. 

Hl = 0.2, ft-fl = 1, the monotonicity parameter H — 0.5, the CFL number Cr = 
Uo^t/'Ax — 0.5, and the parameter xr = 1 in ([37|l . The results of simulation 
are depicted with the exact solution in Figure [21 We note (Figure [2|) that the 
first order scheme, (|37]), exhibits a typical second-order nature, however spuri- 
ous solutions are produced by the scheme. Notice, the numerical simulations 
were performed with such values of the parameter >c, CFL number, Cr, and 
monotonicity parameter, K, that (|57p was not violated. As it can be seen in 
Figure m the boundary maximum principle is not violated by the scheme, i.e., 
the maximum positive values of the dependent variable, v, occur at the bound- 
ary t = 0. It is interesting that the spurious solution (see Figure [2]) produced 
by the scheme COSl has the monotonicity property [19], since no new local ex- 
trema in x are created as well as the value of a local minimum is non-decreasing 
and the value of a local maximum is non-increasing. 

Let us note that the problem of building free-of-spurious-oscillations schemes 
is, in general, unsettled up to the present. Even the best modern high-resolution 
schemes can produce spurious oscillations, and these oscillations are often of 
ENO type (see, e.g., [43] and references therein). We found that the oscillations 
produced by the COSl scheme, ([57]). are of ENO type, namely their amplitude 
decreases rapidly with decreasing the time-increment At, and the oscillations 
virtually disappear under a relatively low CFL number, Cr < 0.15. However, 
the reduction of the CFL number causes some smearing of the solution. The 
spurious oscillations (see Figure [2]) can be eradicated without reduction CFL 
number, but decreasing the parameter >c. Particularly, the spurious oscillations 
disappear if x — 2/3, Cr — 0.5, however, this introduces more numerical 
smearing than in the case of the CFL number reduction. Satisfactory results 
are obtained under >c = 0.82 {Cr — 0.5). The results of simulations are not 
depicted here. 

To gain insight to why the scheme COSl, (|37p . can exhibit spurious solutions, 
let us consider the, so called, first differential approximation of this scheme ( [131 
p. 45], [ini p. 376]; see also 'modified equations' in T3s p. 45], [32], [35]). As 
reported in [13], [49], this heuristic method was originally presented by Hirt 
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(1968) (see [H P- 45]) as well as by Shokin and Yanenko (1968) (see [13 p. 
376]), and has since been widely employed in the development of stable difference 
schemes for PDEs. 

We found that the local truncation error, -0, for the scheme COSl can be 
written in the following form 

{l-x){^xf d^u{x,t) At d'fju) 
^ AAt dx^ 4 dtdx 

o(^^ + {Atf + {Axfy (83) 

By virtue of (|83p . we find the first differential approximation of the scheme 
COSl 

du df{u) At d ( du{x,t) \ 

where B ^ [l - x) [Ax/Atf - A^. The term in right-hand side of (HH) wiU be 
dissipative if 

(l->^)(^)"-^'>0, =^ Cl <!->.. (85) 

Thus, the scheme COSl, ([57)) . is non-dissipative under k — \, and hence can 
produce spurious oscillations. Notice, if >f = 0.82, then we obtain from (|85p 
that Cr < 0.42. Nevertheless, as it is reported above, satisfactory results can 
be obtained under Cr = 0.5 as well. 

So then, the notion of first differential approximation has enabled us to 
understand that the spurious solutions exhibited by the scheme COSl, (j37p . are 
mainly associated with the negative numerical viscosity introduced to obtain the 
scheme of the second order in space, i.e. 0{{Ax)^ + At). Let us consider the 
scheme C0S2, approximating ([T7)l with the accuracy o{{Axy + (Aty). 
Notice, the second order scheme C0S2, p2|) . is nothing more than the scheme 
COSl, (|57|) . with the additional non-negative numerical viscosity. To test the 
scheme C0S2, (|42l) . the inviscid Burgers equation was solved under the initial 
condition (|78p . The numerical solutions were computed under the same values 
of parameters as in the case of the scheme COSl, but C,. = 1. The results of 
simulation are depicted with the exact solution in Figure [3l 

We note (Figure [3]) that the scheme C0S2, (HH), exhibits a typical second- 
order nature without any spurious oscillations. Increasing the value of H (up 
to H = 1) leads to a minor improvement of the numerical solutions, whereas 
decreasing the value of Cr leads to a mild smearing of the solutions. The 
results of simulations are not depicted here. 
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Figure 3: Inviscid Burgers equation. The scheme C0S2 = 1, >f = 1) ver- 
sus the analytical solution. Crosses: numerical solution; Solid line: analytical 
solution and initial data. = 1, H = 0.5, Ax — 0.01. 



5.2 Hyperbolic conservation laws with relaxation 

Let us consider the model system of hyperbolic conservation laws with relaxation 
developed in ^5] : 

- + -(-u^ + a.j=0, (86) 

where 

Qiw, z) = z — m{u — Wo), u ~ w — qaz, (88) 

T denotes the relaxation time of the system, qq, m, a, and uq are constants. The 
Jacobian, A, can be written in the form 

^ ^ I w - + a -9o (w - q^z) | ^^^^ 

The system ([55)) - (l57)) has the following frozen [55' characteristic speeds Ai = a, 
A2 = u + a. The equilibrium equation for ([5S|) - (|57)l is 



dw d f 1 



dt dx 

where 



—u 



1 



aw = 0, (90) 



m 

= ui — (?o2^*, ^* = — Wo) • (91) 

1 -t- mgo 

The equilibrium characteristic speed A* can be written in the form 

A* iw) ^ *^ ' + a. (92) 
1 + mqo 

Pember's rarefaction test problem is to find the solution {w, z} to (l86l) - ((87)) . 
and hence the function u = u{x,t), under t 0, and where 

{w,z} = l {^L,z.in.L)h x<xo (33) 

1 {WR,Z^['Wii,)\ , X>Xo 
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< ul = wl - go^* (wl) < ur = wr - qoz* (wr) . (94) 

The analytical solution of this problem can be found in [15] . The parameters of 
the model system are assumed as follows: qo — —1, m = —1, mq = 3, a = ±1, 
T = 10~^. The initial conditions of the rarefaction problem are defined by 

ul = 2, =^ zl = m {ul - uq) = 1, wl = ul + Qozl = 1, (95) 

ur = 3, =^ zr = m (ur - Uq) = 0, Wr ^ur + qoZR = 3. (96) 

The position of the initial discontinuity, xq, is set according to the value of a 
so that the solutions of all the rarefaction problems are identical [15]. Let a 
position, z^, of leading edge or a position, x\^, of trailing edge of the rarefaction 
be known (e.g., = 0.85, a;^ = 0.7 in f5S]), then 

^0 = x'r - ( +a)t^xi~ ( + a) t. (97) 

At t = 0.3, under dHn])-® we have [15| 

X < 0.7 

0.7 < x < 0.85 . (98) 
X > 0.85 

The results of simulations, based upon the scheme C0S2, ([42]) . together with 
(|76p. under different values of the parameter a(a = l,a — — 1) and different 
values of a grid spacing. Ax, are depicted in Figure |U 

One can clearly see (Figure S]) that the scheme C0S2 is free from spurious 
oscillations. Let us also note that the results generated by the scheme C0S2 are 
less accurate in the case of negative value of a than those in the case of positive 
value of a. Specifically, in the numerical solutions produced under a = — 1, 
the representations of the trailing and leading edges of the rarefaction are more 
smeared than those in the solutions produced under a ~ 1. Notice, under some 
negative value of a, the frozen and the equilibrium characteristic speeds do not 
all have the same sign. 
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0.85- 
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5.3 1-D Euler equation of gas dynamics 

In this subsection we apply the second order scheme C0S2, (1421) . to the Euler 
equations of gamma-law gas: 

^"a?^^ +^F(u) = 0, xGR, t>0; u (x, 0) = u" (x) , (99) 

u ={mi,U2,U3}^ = {p,pt',e}^ , Y [u) = [pv.pv^ +p,{e + p)v]^ , (100) 

P 1 

e = h —pv^^ 7 ~ const, (101) 

7 — 1 2 
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Figure 4: Pember's rarefaction test problem. The second-order scheme C0S2 
= 1, >f = 1) versus the analytical solution for u. Dashed line: numerical 
solution; Solid line: analytical solution. Time t — 0.3, Courant number Cr = 1, 
monotonicity parameter H = 1. (al): Ax = 10^^, a = 1; (a2): Ax = 2.5 x 10^"*, 
a = 1; (bl): Ax = 10"^, a = -1; (b2): Ax = 2.5 x lO""*, a = -1. 



where p, v, p, e denote the density, velocity, pressure, and total energy respec- 
tively. We consider the Riemann problem subject to Riemann initial data 

u°(a;)^("^ a;<2;o ul,ur = const. (102) 

The analytic solution to the Riemann problem can be found in [32l Sec. 14]. 

First we solve the shock tube problem (see, e.g., [5], [32], [33]) with Sod's 
initial data: 









' 0.125 ] 


UL = < 














0.25 J 



(103) 



Following Balaguer and Conde [5] as well as Liu and Tadmor [33] we assume 
that the computational domain is < < 1; the point xq is located at the 
middle of the interval [0, 1], i.e. xq = 0.5; the equations (|99l) are integrated up 
to i = 0.16 on a spatial grid with 200 nodes as in ^ and in [35 • The CFL 
number is taken to be Cr = 1 in contrast to [5] and [33j , where the simulations 
were done under At = O.lAx (i.e. 0.13 < Cr < 0.22). The results of simulations 
are depicted in Figure [5] 

The results depicted in Figure [5] (left column) are not worse in comparison 
to the corresponding third-order central results of [33j p. 418] as well as to 
the results obtained by the fourth-order non-oscillatory scheme in [5] p. 472]. 
Notice, the fourth-order scheme [5] p. 472] gives a better resolution but, in 
contrast to the scheme C0S2, can produce spurious oscillations. 
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Figure 5: Sod's problem. The scheme C0S2 under = 1, H = 0.5, ^ = 1, k = 
0.8 versus the analytical solution. Time t = 0.16, spatial increment Aa; = 0.005 
(left column) and Aa; = 0.0025 (right column). 



24 



Let us also note that the number of multiplications and divisions per one 
grid node in Scheme C0S2, (|42)) . and in the second-order scheme considered 
in [33j is approximately the same, but less than this number in the third- and 
fourth-order schemes developed in [331 P- 418], [SJ p. 472], respectively. Hence, 
the results depicted in Figure [5] (right column) is rather chipper (in terms of 
CPU time) than the ones demonstrated in [331 P- 418], [5,, p. 472]. Along 
with loss of computational efficiency, simulations with low CFL number can, 
in general, lead to excessive numerical smearing. As it is demonstrated above. 
Scheme C0S2 is free from such drawbacks. 



5.4 3-D axial symmetric gas dynamics 

We consider an adiabatic expansion of a gas plume into vacuum [3], i.e., the 
so called Anisimov's problem. Taking into account the symmetry of the plume 
with respect to the axis z, the gas-dynamic equations can be written (for < r, 
z < oo) as follows. 

dp Id (rpvr) d [pvz 



dt r dr dz 



= 0, 



d , . 1 9 r 



dt r dr 

d , . d 



[pVr) + -Tr rp [Vrf + ^ ipVzVr) + „ 

„ ^„ J dr 



d .dp 



dt " d 



z 



19 -.dp 



r 



dr dz 





(104) 


0, 


(105) 


0, 


(106) 


0. 


(107) 




(108) 



dpE Id, , ^ d , , ^ 

-^ + -jr[rvr {pE+p)] + — [v, {pE + p)] = 
dt r dr dz 

p 

pE = h O.dpv^, v'^ ~ v'^ + v^, J — const 

7-1 

The initial conditions are the following (in details, see [3!) 

p — p {r, z) , p/ p'^ = const, Vr — = 0, r, z > 0, t — 0. (109) 

At r = we assume that the axis z is a reflection line. It prohibits any normal 
flux of mass through the boundary r = 0, i.e. 

Vr = 0, r = 0, z > 0. (110) 

Moreover, it is assumed that the pressure (p), density (p), and tangential velocity 
(vz) are even functions of normal distance to the axis z while the normal velocity 
(vr) is an odd function of r. It is also assumed that the plane z = is a reflection 
surface, i.e. the pressure (p), density (p), and tangential velocity (vr) are even 
functions of normal distance above the target surface while the normal velocity 
(vz) is an odd function of z. The analytic solution to the problem p04p - (|110p 
can be found in [3]. 

Notice, every point on the axis r = is a singular point for System (|104p - 
P07p . Assuming that all terms at (|104p are bounded values at a vicinity of 
r = 0, we find that Vr ^ as r ^ 0. Hence 

lim P^'^'^ryO ^ jjj^ P^rUo- Pt'rl.^o ^ d jpVr) ^ 

r-»o+o r r dr ' 
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In perfect analogy we obtain 



p{Vr) dp{Vr)^ pVzVr dpV^Vr Vr {pE + p) dVr {pE + p) 



dr 



dr 



dr 



(112) 



as r ^ 0. By virtue of pTT)) - pT^ . we obtain from (fTM)l - (fTU7|) the following 
conditions at r = 0: 



dp 

dt 



d 



d 



+ — {2pvr) + — {pv,) = 



dr 



dz 



d_ 

dr 
d_ 

dz 



— [pv^Vr) = 0, 

dz 



dpE 



d_ 

dr 



2p{Vr) +p 
P{Vzf +P 

[2vr {pE +p)] + ^ [vz {pE + p)] = 0. 



— {2pVzVr) = 0, 

dr 



d 



(113) 

(114) 
(115) 
(116) 



dt dr dz 
In the analytic solution [5] of the problem (|104p - (|110p the following input 
data are required: the initial dimensions of the plume, Rq and Zq, its mass Mp, 
and the initial energy Ep. We will use the following values as the reference 



Ro, = v/(57 - 3) Ep/Mp, = U/v^, p^ 



Mf 



{RnZo) 



quantities: U 

p* = p^vl. 

The equations p04p - ()107p are integrated up to t = 0.4 with a = Zq/^Rq = 
0.1. The CFL number is taken to be Cr — 1. It is assumed that the spatial 
increments are the following: Ar = 0.0025, Az = 0.00025 if < t < 0.05; 
Ar = 0.0025, Az = 0.0005 if 0.05 < t < 0.1; Ar = 0.005, Az = 0.001 if 
0.1 < i < 0.2; Ar = 0.01, Az = 0.002 if 0.2 < i < 0.4. The results of 
simulations as well as the analytical solution are depicted in Figures [H [71 We 
observe (Figures [6l [7]) that the numerical and analytical solutions are practically 
coincide, but in the vicinity of the front, namely, for very small values of density. 

6 Appendix 

Proposition 5 Let us find the order of accuracy, r, in {21^ if di will be ap- 
proximated by di with the order of accuracy s, i.e. let 



d, = d, + 0{{Axf). 
Let U {x) be sufficiently smooth, then we can write 



(117) 

(118) 
(119) 



26 




Figure 6: Anisimov's problem, density and pressure distribution. C0S2 scheme 
versus the analytical solution, a = Zq/Rq = 0.1, time t = 0.4, CFL number 
Cr = 1, monotonicity parameter H = x = 1, spatial increments: Ar = 0.0025, 
Az = 0.00025 if < f < 0.05; Ar = 0.0025, Az = 0.0005 if 0.05 < t < 0.1; 
Ar = 0.005, Az = 0.001 if O.K t < 0.2; Ar = 0.01, Az = 0.002 if 0.2 < t < 0.4. 
Dashed lines: numerical solution; Solid lines: analytical solution. 
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N0.3-: 





N0.3- 



0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 



Figure 7: Anisimov's problem, momenta {pVz and pVr) distribution. C0S2 
scheme versus the analytical solution, a = Zq/Rq = 0.1, time t = 0.4, CFL 
number Cr = 1, monotonicity parameter W = >c = 1, spatial increments: Ar = 
0.0025, Az = 0.00025 if < f < 0.05; Ar = 0.0025, Az = 0.0005 if 0.05 < 
t < 0.1; Ar = 0.005, Az = 0.001 if 0.1 < i < 0.2; Ar = 0.01, Az = 0.002 if 
0.2 <t< 0.4. Dashed lines: numerical solution; Solid lines: analytical solution. 
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Combining the equalities I1118\) and \119[ we obt 



am 



Q2JJ n2 



In a similar manner we write: 



i+05 



+0((Axf). (120) 



2 



= + c/;;o5— + ( — ) + o ( {^xf ) , (121) 



rf. = ^.:+05 - c^;;o5^ + \u'^u, (^)' + ^ ((ax)=^) . (122) 



Subtracting the equations hl21]) and il22\) . we obtain 



i+05 



o 



[(Axf) . (123) 



In view of il23\) and |j_?7p we obtain from il20\) the following interpolation 
formula 

C/.+05 - I (f/.+i + - di) + O ((Ax)^ + (Ax)^+') . (124) 

In view of \124^ we obtain that r = min (4, s + 1) . 

Proposition 6 Let us construct a second order scheme based on operator- 
splitting techniques. We will, in fact, use the summarized (summed) approx- 
imation method J4^\ Section 9.3] to estimate order of approximation. Consider 
the following equation 

Vu = Viu + T'su = ^ - Lu - 0, VkU=--^~ LfcU, k - 1, 2, (125) 

where Lk is an operator, e.g. a differential operator, a real analytic function, 
etc., acting on u [x, t). We approximate I1125\) on the cell [xi-i, Xi^i] x [t„, 
by the following difference equation with the accuracy 0{{Ax) + (At) ) 

v"+l - v" 

nv = ' ' - Aiv"+"-5 - A2v"+°-5 = 0, (126) 

where it is assumed that the operator L^u is approximated by the operator A^u 
with the accuracy 0{(Ax)'^), i.e. 

A,u"+°-'^ = (L,u)f+°-^ + O [{Axf) . (127) 

In view of the operator splitting idea, to the problem hl26\) there corresponds the 
following chain of difference schemes 

1 w"+"-^ - w" 

"^-"2 '0.5A. ' -Aiwr-^-O, (128) 
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"2 O.SAt ^^^"2 v^--; 

One can see from the above that the operator VkU is approximated by IlfeU with 
the accuracy 0{At + (Aa;)^) 

ni<+«-^ = (Piu)r°-^ - ^ (^1^ j + O {{Atf + {Axf) , (130) 

n^u^o-^ = {V,u)r'-' + ^ (^) + O {{Atf + (Axf) . (131) 

In view of U30\} - !!73l]) . the local truncation error \32\ p. 142], ip, on a suffi- 
ciently smooth solution u{x, t) to \125]) is found to be 

V' = nu = Hiu+nzu = (132) 

(Piu+P2u)"+" ' + O ((At)' + (A:e)') = O ((At)' + (Ax)') . (133) 

Thus, the implicit scheme, il28]) . together with the explicit scheme, U29\) . ap- 
proximate il25\) with the second order. 
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